#include "field.hpp"
#include "meshTypeFunctions.hpp"
#include "bc_weir_funcs.hpp"
#include "bc_weir_bottomElevation.hpp"
#include "bc_weir_velocity.hpp"
#include "bc_weir_waterdepth.hpp"
#include "util.hpp"

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

using namespace hms;
using Type = hms::boundaryCondition::Type;
using Catch::Approx;

TEST_CASE("weir-like boundary functions"){
	using namespace hms::bc_weir;


	REQUIRE( degToRad(90) == Approx(0.5*math::pi) );

	REQUIRE( pressureFactor(0) == Approx(1) );
	REQUIRE( pressureFactor( degToRad(90) ) == Approx(0.1) );
	REQUIRE( pressureFactor( degToRad(15)  ) == Approx(0.933012701892219) );
	REQUIRE( pressureFactor( degToRad(75) ) == Approx(0.1) );

	scalar alpha, cp;
	Array1Ns<1> d_e, d_e_target, d, v;

	auto test = [&](scalar targetValue){
		waterDepth_edge(d_e,d,v,cp);
		scalar d_e_val { d_e[0] };
		REQUIRE( d_e_val == Approx(targetValue) );
	};

	d = 0.2;
	v = 1;
	alpha = degToRad(15);
	cp = pressureFactor(alpha);
	test(0.121863060849659);

	v = 2;
	test(0.194316422803081);

	alpha = degToRad(75);
	cp = pressureFactor(alpha);
	test(0.162732159227197);

	v = 1;
	test(0.102283899515963);
}

TEST_CASE("weir-like boundary conditions"){
	auto test = [&](auto&& dummyMesh){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;

		auto mesh { uniMeshAs<MeshType>(1,1,1,1,1,true) };
		const BoundaryPatch
			& bottom { mesh.boundaryPatch("bottom") },
			& right  { mesh.boundaryPatch("right" ) };
		
		Field
			v { mesh, 2 },
			d { mesh, 1 },
			z { mesh, 1 };

		d.values() = 0.2;
		z.values() = 0;
		v.values() = 1;

		auto setBottomVel = [&](scalar vOrth, scalar vPara){
			v.domainValues(bottom) = Array2s{ vPara, -vOrth };
		};
		auto setRightVel = [&](scalar vOrth, scalar vPara){
			v.domainValues(right) = Array2s{ vOrth, vPara };
		};

		SECTION("water depth"){
			/* bottom boundary, requires coord transformation of velocity */
			setBottomVel(1,2);

			auto& bottomBC { d.setBoundaryCondition(
				Weir_waterdepth{d, bottom, v, degToRad(15)}
			) };
			REQUIRE( bottomBC.alpha() == degToRad(15) );
			REQUIRE( bottomBC.cp() == bc_weir::pressureFactor( degToRad(15) ) );
			bottomBC.update();
			REQUIRE( d.ghostValues(bottom)(0,0) == Approx(0.1) );

			bottomBC.alpha( degToRad(75) );
			bottomBC.update();
			REQUIRE( d.ghostValues(bottom)(0,0) == Approx(0.1) );

			setBottomVel(2,1);

			bottomBC.update();
			REQUIRE( d.ghostValues(bottom)(0,0) == Approx(0.125464318454394) );

			bottomBC.alpha( degToRad(15) );
			bottomBC.update();
			REQUIRE( d.ghostValues(bottom)(0,0) == Approx(0.188632845606162) );

			/* right boundary, no transformation required */
			setRightVel(1,2);

			auto& rightBC { d.setBoundaryCondition(
				Weir_waterdepth{d, right, v, degToRad(15)}
			) };
			rightBC.update();
			REQUIRE( d.ghostValues(right)(0,0) == Approx(0.1) );

			rightBC.alpha( degToRad(75) );
			rightBC.update();
			REQUIRE( d.ghostValues(right)(0,0) == Approx(0.1) );

			setRightVel(2,1);

			rightBC.update();
			REQUIRE( d.ghostValues(right)(0,0) == Approx(0.125464318454394) );

			rightBC.alpha( degToRad(15) );
			rightBC.update();
			REQUIRE( d.ghostValues(right)(0,0) == Approx(0.188632845606162) );
		}
		SECTION("velocity"){
			/* bottom boundary, requires coord transformation of velocity */
			setBottomVel(1,2);

			auto& bottomBC { v.setBoundaryCondition(
				Weir_velocity{v, bottom, d, degToRad(15)}
			) };
			bottomBC.update();
			REQUIRE( v.ghostValues(bottom)(0,0) == Approx(2) );
			REQUIRE( v.ghostValues(bottom)(1,0) == Approx(-5./3) );

			bottomBC.alpha( degToRad(75) );
			bottomBC.update();
			REQUIRE( v.ghostValues(bottom)(0,0) == Approx(2) );
			REQUIRE( v.ghostValues(bottom)(1,0) == Approx(-5./3) );

			setBottomVel(2,1);

			bottomBC.update();
			REQUIRE( v.ghostValues(bottom)(0,0) == Approx(1) );
			REQUIRE( v.ghostValues(bottom)(1,0) == Approx(-2.91605349427636) );

			bottomBC.alpha( degToRad(15) );
			bottomBC.update();
			REQUIRE( v.ghostValues(bottom)(0,0) == Approx(1) );
			REQUIRE( v.ghostValues(bottom)(1,0) == Approx(-2.11699633237235) );

			/* right boundary, no transformation required */
			setRightVel(1,2);

			auto& rightBC { v.setBoundaryCondition(
				Weir_velocity{v, right, d, degToRad(15)}
			) };
			rightBC.update();
			REQUIRE( v.ghostValues(right)(0,0) == Approx(5./3) );
			REQUIRE( v.ghostValues(right)(1,0) == Approx(2) );

			rightBC.alpha( degToRad(75) );
			rightBC.update();
			REQUIRE( v.ghostValues(right)(0,0) == Approx(5./3) );
			REQUIRE( v.ghostValues(right)(1,0) == Approx(2) );

			setRightVel(2,1);

			rightBC.update();
			REQUIRE( v.ghostValues(right)(0,0) == Approx(2.91605349427636) );
			REQUIRE( v.ghostValues(right)(1,0) == Approx(1) );

			rightBC.alpha( degToRad(15) );
			rightBC.update();
			REQUIRE( v.ghostValues(right)(0,0) == Approx(2.11699633237235) );
			REQUIRE( v.ghostValues(right)(1,0) == Approx(1) );
		}

		SECTION("bottom elevation"){
			/* bottom boundary, requires coord transformation of velocity */
			setBottomVel(1,2);

			auto& bottomBC { z.setBoundaryCondition(
				Weir_bottomElevation{z, bottom, d, v, degToRad(15)}
			) };
			bottomBC.update();
			REQUIRE( z.ghostValues(bottom)(0,0) == Approx(0.0223338672879957) );

			bottomBC.alpha( degToRad(75) );
			bottomBC.update();
			REQUIRE( z.ghostValues(bottom)(0,0) == Approx(0.0223338672879957) );

			setBottomVel(2,1);

			bottomBC.update();
			REQUIRE( z.ghostValues(bottom)(0,0) == Approx(-0.131411453193722) );

			bottomBC.alpha( degToRad(15) );
			bottomBC.update();
			REQUIRE( z.ghostValues(bottom)(0,0) == Approx(-0.0128291122563493) );

			/* right boundary, no transformation required */
			setRightVel(1,2);

			auto& rightBC { z.setBoundaryCondition(
				Weir_bottomElevation{z, right, d, v, degToRad(15)}
			) };
			rightBC.update();
			REQUIRE( z.ghostValues(right)(0,0) == Approx(0.0223338672879957) );

			rightBC.alpha( degToRad(75) );
			rightBC.update();
			REQUIRE( z.ghostValues(right)(0,0) == Approx(0.0223338672879957) );

			setRightVel(2,1);

			rightBC.update();
			REQUIRE( z.ghostValues(right)(0,0) == Approx(-0.131411453193722) );

			rightBC.alpha( degToRad(15) );
			rightBC.update();
			REQUIRE( z.ghostValues(right)(0,0) == Approx(-0.0128291122563493) );
		}
	};

	SECTION("StructMesh"){
		test(StructMesh{});
	}
	SECTION("RectMesh"){
		test(RectMesh{});
	}
	SECTION("UniMesh"){
		test(UniMesh{});
	}
	SECTION("Mesh"){
		test(Mesh{});
	}
}
